Deciphering Core Collapse Supernovae: Is Convection the Key? 
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ABSTRACT 

We couple two-dimensional hydrodynamics to detailed one-dimensional multigroup 
flux-limited diffusion neutrino transport to investigate prompt convection in core 
collapse supernovae. Our initial conditions, time-dependent boundary conditions, and 
neutrino distributions for computing neutrino heating, cooling, and deleptonization 
rates are obtained from one-dimensional simulations that implement multigroup 
flux- limited diffusion neutrino transport and one-dimensional hydrodynamics. The 
development and evolution of prompt convection and its ramifications for the shock 
dynamics are investigated for both 15 and 25 M models, representative of the two 
classes of stars with compact and extended iron cores, respectively. In the absence of 
neutrino transport, prompt convection develops and dissipates on a time scale ~15 
ms for both models. Prompt convection seeds convection behind the shock, which 
causes distortions in the shock's sphericity, but on the average, the shock radius is 
not boosted significantly. In the presence of neutrino transport, prompt convection 
velocities are too small relative to bulk inflow velocities to result in any significant 
convective transport of entropy and leptons. A simple analytical model supports our 
numerical results, indicating that the inclusion of transport reduces the convection 
growth rates and asymptotic velocities by factors of 4-250. 

1. Introduction 

Ascertaining the core collapse supernova explosion mechanism is a long-standing problem in 
astrophysics. The current paradigm begins with the collapse of a massive star's iron core and the 
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generation of an outwardly propagating shock wave that results from core rebound. Because of 
dissociation and neutrino losses the shock stagnates. This sets the stage for a shock reheating 
mechanism whereby neutrino energy deposition via electron neutrino and antineutrino absorption 
on nucleons behind the shock reenergizes it (Bethe & Wilson 1985; Wilson 1985); however, no 
recent numerical simulations of shock reheating and the subsequent shock evolution produce 
explosions unless the neutrino luminosities or the postshock neutrino heating efficiencies are 
boosted by some other mechanism. One mechanism may be convection (Wilson & Mayle 1993, 
Herant et al. 1994, Burrows et al. 1995, Janka & Muller 1996). 

There are at least three modes of convection that may develop during the shock reheating 
phase. (1) Prompt convection near and below the "neutrinospheres," which is the subject of this 
Letter and which may occur immediately after the formation of the shock (Burrows & Fryxell 
1993, Bruenn & Mezzacappa 1994, Janka & Muller 1996). (2) Doubly diffusive instabilities below 
the neutrinospheres [e.g., "neutron fingers" (Wilson & Mayle 1993); however, see Bruenn Sz Dineva 
1996]. (3) Neutrino-driven convection below the shock (Herant et al. 1992, Miller et al. 1993, 
Herant et al. 1994, Burrows et al. 1995, Janka & Muller 1996, Mezzacappa et al. 1996, Calder et 
al. 1996b). 

Prompt convection may initially be driven by negative entropy gradients imprinted on the 
core by the weakening shock and may shortly thereafter be sustained by the negative lepton 
gradient that arises near and below the electron neutrinosphere as a result of electron neutrino 
losses. Based on multidimensional hydrodynamics and simplified neutrino transport, it has been 
suggested that prompt convection increases the neutrinosphere luminosities (Burrows & Fryxell 
1993, Herant et al. 1994, Janka & Muller 1996) and pushes the shock out farther in radius (Janka 
k, Muller 1996). Based on one-dimensional hydrodynamics, mixing-length convection, and detailed 
multigroup flux-limited diffusion (MGFLD) neutrino transport, all coupled self-consistently, it 
has been shown that prompt convection does not lead to a significant increase in the neutrino 
luminosities and has little effect on the shock dynamics (Bruenn & Mezzacappa 1994, Bruenn et 
al. 1995). In this Letter, we remove the restriction to spherical symmetry in the hydrodynamics 
and the need for mixing-length convection, and couple two-dimensional hydrodynamics to detailed 
one-dimensional MGFLD neutrino transport. 



2. Initial Models, Codes, and Methodology 

We began with the 25 M precollapse model S25s7b provided by Woosley (1995). The 
initial model was evolved through core collapse and bounce using MGFLD neutrino transport and 
Lagrangian hydrodynamics. The one-dimensional data at 8 ms after bounce (297 ms after the 
initiation of core collapse) were mapped onto our two-dimensional Eulerian grid. The inner and 
outer boundaries of our grid were chosen to be at radii of 25 km and 1000 km, respectively. 128 
nonuniform radial spatial zones were used. This gave sufficient resolution at the inner boundary 
and at the shock. The nonuniform zoning was varied until no significant differences were seen 
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between one-dimensional runs using 128 nonuniform and 512 uniform zones. 128 uniform angular 
zones spanning a range of 180 degrees and reflecting boundary conditions were used for 0. 

The inner boundary was chosen to be below the unstable region at the onset of the simulation, 
which at 8 ms extended from 33 to 58 km. At the base of the unstable region, s = 6.2 and 
Y e = 0.23, and at the top, s = 4.7 and Y e = 0.16, where s is the entropy per baryon in units 
of Boltzmann's constant and Y e is the electron fraction. Spherically symmetric time-dependent 
boundary data for the two-dimensional hydrodynamics were supplied by our one-dimensional run. 
At each time step, the boundary data at our fixed inner Eulerian radius were extracted from the 
one-dimensional data by interpolation in r and t. The outer boundary data were specified in the 
same way. 

The two-dimensional hydrodynamics was evolved using an extended version, EVH-1, of the 
PPM hydrodynamics code VH-1 developed by Blondin and colleagues at the Virginia Institute for 
Theoretical Astrophysics. Our extension allows for coupling to general equations of state. The 
matter in our simulations was in nuclear statistical equilibrium; to describe its thermodynamic 
state, we used the equation of state provided by Lattimer and Swesty (Lattimer &; Swesty 1991). 

Because the finite differencing in our PPM scheme is nearly noise free and because we cannot 
rely on machine roundoff to seed convection in a time short compared with the hydrodynamics 
time scales in our runs, we seeded convection in the Ledoux unstable regions below and around the 
neutrinospheres by applying random velocity perturbations to the radial and angular velocities 
between ±5% of the local sound speed. 

In our two-dimensional simulations, gravity was assumed to be spherically symmetric. The 
gravitational field in the convectively unstable region was dominated by the enclosed mass at the 
region's base, which at the start of our simulations was 0.82 M . The enclosed mass at the top of 
the unstable region at this time was 1.1 M . Moreover, at t = 303 ms in our simulation without 
neutrino transport (6 ms after the initiation of our run, at a time when prompt convection was 
fully developed), the density varied with 6 about its average value (which was used in computing 
the spherically symmetric gravitational field) between -17% and +24% at 41 km and -9% and 
+57% at 50 km. Therefore, we do not believe this assumption was a serious shortcoming of our 
calculations. The time dependence of the mass enclosed by our inner boundary, given by our 
one-dimensional MGFLD runs, was taken into account. The solution of the Poisson equation for 
the gravitational potential will be incorporated in future investigations. 

The neutrino heating and cooling, and the change in the electron fraction, were computed 
using the following formulae: de/dt = cY,j=i J E^dE u [tpf / \\ a ^ — — Y>°)]/ p(hc) 3 and 
dY e /dt = cniB Ef=i OLi J EldE^JX^ - - i>f)]/p(hcf, where e is the internal energy per 
gram; E v , ip®, \\ a \ and ji are the electron neutrino or antineutrino energy, zeroth distribution 
function moment, absorption mean free path, and emissivity, respectively; m-Q is the baryon mass; 
i = l,cxi = 1 corresponds to electron neutrinos, and i = 2,ct2 = —1 corresponds to electron 
antineutrinos. The time-dependent ijj^s are obtained from tables in r and t constructed from our 



-4 - 



one-dimensional MGFLD simulations. Comparisons were made to ensure that in one dimension the 
results from our code matched the results obtained with Bruenn's MGFLD code, modulo EVH-l's 
better resolution of the shock. 

In the 15 M case, we began with the precollapse model S15s7b provided by Woosley (1995). 
The one-dimensional data at 16 ms after bounce were mapped onto our two-dimensional Eulerian 
grid, and the inner and outer boundaries were chosen to be at radii of 20 km and 1000 km, 
respectively. Initially, the unstable region extended from 29 to 58 km, with s and Y c varying 
from 6.1 to 5.3 and 0.26 to 0.15, respectively. At this time, the enclosed mass at the bottom 
and top of the unstable region was 0.77 M and 1.1 M , respectively. At t = 218 ms in our 
simulation without neutrino transport (10 ms after the initiation of our run, at a time when 
prompt convection was fully developed), the angular density variations about the average density 
were between -15% and +18% at 41 km and -13% and +65% at 51 km. 

3. Results 

For both our 15 and 25 M models, two simulations were carried out for a duration of 100 
ms after the initial postbounce time — one with neutrino transport and one without. 

For our 25 M model without transport, prompt convection develops and dissipates in ~15 
ms, by which time the initial gradients are smoothed out by convection. The entropy evolution is 
shown in Figure 1. Prolonged convection would require the maintenance of the electron fraction 
gradient by electron neutrino escape (neutrino transport) at the neutrinosphere. For the duration 
of the prompt convection episode, the fluid velocity at the inner boundary remained positive, 
ensuring that none of this innermost convection was swept off our grid. 

Prompt convection seeds unstable regions at successively larger radii and eventually 
convection reaches the shock and distorts it from sphericity. However, on the average, there is no 
significant difference between our one- and two-dimensional shock trajectories: convection does 
not significantly boost the stalled shock radius. [For details, see Calder et al. (1996b).] 

The outward motion of the shock in Figure 1 is not indicative of an explosion but results from 
the decreasing preshock accretion ram pressure and is seen in our one-dimensional run. 

With transport, both the convection growth rate and the asymptotic convection velocities 
are substantially smaller (see Section 4). The asymptotic velocities become too small relative to 
the bulk inflow to result in any significant convective transport of entropy and leptons, and the 
evolution proceeds as it does in our one-dimensional run. The evolution is shown in Figure 2. 

The results for the 15 M model are essentially the same. In the absence of neutrino 
transport, convection develops and dissipates in ~15 ms. Prompt convection seeds successive 
unstable regions between the prompt convection region and the shock. The shock is distorted from 
sphericity, but again, on the average, is not boosted significantly in radius relative to its position 
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in our one-dimensional run without convection. 

As in the 25 M case, the prompt convection growth rate and asymptotic velocities in the 
presence of neutrino transport are too small to result in any significant convective transport of 
entropy and leptons. Complete details for both models will be given in Calder et al. (1996a). 

For our 15 M model, the velocity at the inner boundary was negative during the course of the 
prompt convection episode; however, this did not have a significant effect on the development of 
convection in our "hydrodynamics only" run. Moreover, the key outcome, that neutrino transport 
inhibits the development of prompt convection, remained unchanged. 



4. Analytical Model 

Convection near or below the neutrinospheres can be profoundly influenced by the neutrino 
transport of energy and leptons between a convecting fluid element and the background. In effect, 
convection becomes "leaky," and differences between a convecting fluid element's entropy and 
lepton fraction and the background's entropy and lepton fraction, from which the buoyancy force 
driving convection arises, are reduced. To construct the simplest model of this, we will assume 
that the lepton fraction gradient is zero and that convection is driven by a negative entropy 
gradient that is constant in space and time. [Reversing the roles of the entropy and lepton 
fraction gradients would give analogous results. The more general case in which both gradients 
are nonzero has been considered by Bruenn and Dineva (1996). This complicates the analysis and 
can lead to additional modes of instability, such as semiconvection and neutron fingers, which are 
not relevant here.] We will also assume that the effect of neutrino transport is to equilibrate the 
entropy of a fluid element with the background entropy in a characteristic time scale t s . If the 
fluid element and the background are in pressure balance, and if we neglect viscosity, the fluid 
element's equations of motion are 

v = -a s 6 s (1) 

P 

9s = -V s - Tr V (2) 

where 9 S = s — s, with s and s being the fluid element's entropy and the background's 
entropy, respectively; g is the local acceleration of gravity; v is the radial velocity; and 
a s = — (d In p/d In s) PYl > 0, where Yi is the common lepton fraction. Equation (|l|) equates the 
fluid element's acceleration to the buoyancy force arising from the difference between its entropy 
and the background's entropy; equation (|2|) equates 6 S to s minus s, where s results from the fluid 
element's equilibration with the background, and s results from its motion through the gradient 
in s. 
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If we neglect neutrino effects (r s = oo), the solutions to equations (|lj) and (Q) indicate that (a) if 
ds/dr > 0, the fluid element oscillates with the Brunt- Vaisala frequency u>bv = [{got s ds /dr) / p] 1 ^ 2 , 
and (b) if ds/dr < 0, it convects, i.e., its velocity increases exponentially, and the convection 

1 /2 

growth time scale is given by r = tbv = [— (gotsds /dr)/p] ■ When neutrino transport effects 
are included in the convectively unstable case (ds/dr < 0), the fluid element convects, but the 
convection growth time scale r > tbv is given by 1/r = + l/^r 2 ] 1 ^ 2 — 1/2t s . In the limit 

t s "C tbv, the growth time scale increases by tbv/ t s, i-e., r ~ t bv /t s . 

In addition to reducing convection's growth rate, neutrino transport also reduces its 
asymptotic velocities. In particular, in the limit t s <C tbv-, the solutions to equations ([[]) and (||) 
show that a fluid element's velocity after moving a distance I from rest is reduced by the factor 
t s /tbv- 

To apply this analysis to prompt convection, we note that tbv is 2 - 3 ms in the region 
between 10 11 and 10 12 g cm -3 in models S15s7b and S25s7b immediately after shock propagation. 
On the other hand, our e from neutrino heating and cooling gives values for r s that decrease from 
0.6 ms at 10 11 g cm -3 to 0.01 ms at 10 12 g cm~ 3 . Our Y e gives values for the lepton equilibration 
time scales that are about a factor of 5 smaller. These imply that neutrino transport should 
reduce the growth rate and asymptotic velocities of entropy-driven convection by factors of from 
~ 4 near 10 11 g cm -3 to ~ 250 near 10 12 g cm -3 . The lepton-driven convection growth rate and 
asymptotic velocities should be reduced by an additional factor of 5. 



5. Caveat 

For model S15s7b, the Planck-averaged optical depth from the top to the bottom of the 
unstable region varies from 1.2 to 24 for electron neutrinos and from 0.56 to 11.3 for electron 
antineutrinos. The corresponding quantities for S25s7b vary from 3.0 to 13.3 and from 1.4 to 5.7, 
respectively. 

Because we are imposing a background neutrino distribution in our two-dimensional 
simulations, in optically thick regions we are overestimating the rate for a fluid element to 
equilibrate with the background, and therefore, overestimating the effect transport has on 
inhibiting the development of prompt convection. The equilibration would be affected by the 
advected trapped neutrinos and by the finite time for neutrino transport between the element and 
its surroundings, both of which are neglected in our analysis. 

Equilibration experiments (Bruenn et al. 1995, Bruenn & Dineva 1996), which include these 
effects, show that a fluid element of one pressure scale height in radius will equilibrate in entropy 
with a time scale of 0.45 ms at 3 x 10 11 g cm~ 3 and 1 ms at 1 x 10 12 g cm -3 . Smaller modes will 
equilibrate faster. The equilibration time scale for Y e is 0.1 - 0.4 ms at 3 x 10 11 g cm -3 and 0.2 - 1 
ms at 1 x 10 12 g cm -3 . These time scales are small compared with the Brunt- Vaisala time scales 
- implying that our conclusions are valid — but not as small as our heating and cooling rates 
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6. Summary 

Near and below the neutrinosphere, neutrino transport equilibrates a convecting fluid element 
with its surroundings in both entropy and electron fraction in a fraction of a millisecond. As 
a result, prompt convection growth rates and asymptotic velocities are reduced by factors of 
4-250. Prompt convection velocities become too small relative to the bulk inflow to result in any 
significant convective transport of entropy and leptons to the neutrinospheres; therefore, prompt 
convection will have no effect on boosting the neutrinosphere luminosities nor on boosting the 
neutrino reheating of the stalled supernova shock wave. 
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Fig. 1. — The 25 M model: without neutrino transport, the evolution in entropy shows the 
development of prompt convection and the convection it seeds between the prompt convection 
region and the shock. 

Fig. 2. — The 25 M model: with neutrino transport, the evolution in entropy shows no significant 
convection. 



